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Direct demodulation method (DDM) was applied to reconstruct y-ray spectra. Boosted Richardson-Lucy 
iteration was introduced into DDM. Monte Carlo method (here GEANT 4) was proposed to calibrate response 
function and establish response matrix. First, gauss function was regarded as total energy peak. Spectra line 
was simulated with nine gauss functions. And afterwards DDM was applied to reconstruct the simulated spectra 
line and determine peak positions and areas. Compared with original spectra, for case that peak position interval 
was about 1/3 full width half maximum (FWHM), the error of rebuilding peak position was 2 channels. The rest 
of peaks could be searched accurately. The relative errors of all peaks’ area were less than 4%. Then, three key 
factors, including noise, background, response matrix, were discussed. Finally, DDM was applied to calibrate 
the field Nal gamma spectrometer. The errors of U, Th, K were less than 5%. Comprehensive studies have 
shown that it is feasible to reconstruct gamma-ray spectra with DDM. DDM can significantly pseudo-improve 
energy resolution of gamma spectrometer, effectively decompose doublets whose peak potential interval is 
1/3 FHWM, and accurately search peak and calculate areas. DDM can restrain noise strongly but is greatly 
influenced by background. And DDM can improve the accuracy of qualitative and quantitative analysis in 


combination with the conventional spectrum analysis method. 
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I. INTRODUCTION 


Gamma energy spectrum analysis is a technique of us- 
ing gamma spectrometer to measure y-rays of radionuclides 
and determine directly categories and contents of particu- 
lar radionuclides by processing the measured y-ray spec- 
tra. Generally, the data processing includes smoothing, back- 
ground deduction, peak searching, calculation of peak ar- 
eas and energy calibration. Certainly, peak position and 
peak area are the key to qualitative analysis and quantita- 
tive analysis respectively. Traditional peak searching meth- 
ods [1] includes IF function method, Gaussian product func- 
tion method, derivative method, covariance method, and sym- 
metrical zero area transformation method. Although they 
can be of high precision for strong radioactive y-ray spec- 
tra and singlet, these peak searching methods are rough for 
weak radioactive y-ray spectra and doublets. Especially for 
weak radioactive y-ray spectra, IF function method, Gaussian 
product function method and covariance are not desirable [2], 
while derivative and symmetric zero area methods have lim- 
ited effects. 

Most of the traditional peak area methods, such as total 
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peak area (TPA) method, Covell peak area method, Was- 
son peak area method, Sterlinski peak area method, Wasson- 
Sterlinski peak area method, Quittner peak area method, W- 
S-Q peak area method, and Q-S peak area method, request de- 
termination of the peak position and border, and can be only 
applied for singlet only [3]. 


The curve fitting method performs better [4, 5]. Its basic 
principle is as follows. A mathematical model, or a func- 
tion, is established to describe the shape of total energy peak 
and baseline, then parameters of the function are determined 
from the experimental data, and finally the peak area is ob- 
tained by integrating the function. In practice, however, due 
to complexity of detectors and detection environment, mea- 
sured y-ray spectra are often too complicated to determine 
the function. In other words, application of the curve fitting 
method of peak area is extremely circumscribed. 


In 1997, de-convolution [6] was proposed to process y-ray 
spectra. For ill-posed problems, conventional de-convolution 
methods, such as maximum entropy method [7], is sensitive 
to the error of input data, namely, a small error can cause 
great oscillation. Fortunately, Tikhonov et al. studied ill- 
posed problems in 1977 on a strict mathematic basis by in- 
troducing the regularization theory and methods [8]. The 
idea is to transform ill-conditioned problems into well-posed 
problems, so as to ensure the solution acceptable and sta- 
ble enough. Regularization falls into three categories gen- 
erally: the least squares, smoothing methods, and iterative. 
The least-squares method is vulnerable to premature and the 
solutions have large error [9]. The degree of smoothing is 
not easy to control and smoothing itself has error [10]. In the 
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iterative approach, one attempts to generate successive ap- 
proximations which converge to an appropriate solution. Per- 
haps the most significant feature of iterative method, for the 
unfolding application, lies in the ability to incorporate read- 
ily most of the major physical implication of the description. 
Also, the necessity for an explicit calculation of the inverse 
matrix is avoided [11]. This is important, as the most response 
matrices are ill-conditioned, and the error in the elements of 
the inverse matrix can be prohibitive. 

There are three commonly used regularization (iteration) 
de-convolution methods: Gold [12], Richardson-Lucy [13], 
and the maximum posteriori algorithms [14]. The direct de- 
modulation method (DDM) [15, 16] was proposed by Li Tibei 
and Wu Mei. Its basic principle is to use some known phys- 
ical conditions to control iterative process for solving mod- 
ulation equation. The Richardson-Lucy iteration converges 
to the maximum likelihood solution. Compared with con- 
ventional de-convolution methods, DDM can be rarely re- 
stricted by shape of y-ray spectra, radioactive intensity, and 
error of input, by taking full advantage of known information 
and utilizing nonlinear physical constraint conditions. DDM 
has succeeded in high-energy astronomy for reconstructing 
astrophysical images of low signal-to-noise ratio (SNR), low 
statistics and low-resolution [17—19]. But there are rarely re- 
ports about DDM for y-ray spectra. In this paper, we propose 
to utilize DDM to reconstruct y-ray spectra. Monte Carlo 
simulation is performed to calibrate response functions and 
establish response matrix. The purpose of this study is to 
improve the energy resolution of gamma spectrometry, and 
accuracy of y-ray spectra analysis. 


Il. THE DIRECT DEMODULATION METHOD 


The input and output of gamma spectrometer can be re- 
garded as a complete system. The incident radiation can be 
imaged as a sum of ó-functions. So at the system output 
the spectrum represents a linear combination of ó-functions 
with various amplitudes located at various channels function, 
which are blurred by the system response functions. The pur- 
pose of DDM was to eliminate the influence of response func- 
tion to the utmost. Ideally, a complete y-ray spectra consist- 
ing of ó-functions can be obtained. 

Let us suppose that the intensity distribution function of 
source is f(a). From a signal processing point of view, mea- 
surement process of spectrum can be regarded as the modula- 
tion process of spectrum of radioactive source. The measured 
y-ray y(j) is the modulated result. The modulation equation 
can be expressed as 


I A(,2)f (z)dz = y(3), (b 


where h(j, x) represents the system response function. 
The matrix equation can be expressed as 


HF =Y, (2) 


where H is system response matrix. For y-ray spectra, H is 
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the lower triangular Toeplitz matrix: 


PE Dono] M 
h(N) h(N -A 


Equation (2) means that spectrum of radioactive source F 
is modulated by gamma spectrometer, and the measured spec- 
trum Y is obtained. The reconstruction process is an inverse 
process: demodulate the spectrum Y to reconstruct spectrum 
of radioactive source F. 

For gamma spectrometry, y(j) and f(x) are both nonnega- 
tive discrete points. It should be emphasized that this is vitally 
important for DDM. Besides, h(j, x) is discretized in the nu- 
merical calculation. Hence the discrete form of Eq. (1) can be 
expressed as 


XO h( A) FG) = ylk), (4) 
k 


where h(k, j) is the element at the k'' row and j column of 
response matrix H, f(j) is the j" element of F and y(k) is 
the k element of Y. However, Eq. (4) is an ill-posed prob- 
lem. So DDM is applied for Eq. (4). 

Richardson-Lucy iteration algorithm (R-L) is based on 
Bayesian statistical theory. It can be expressed as 


2, Ms 


where d? (k) = Y> h(k, j) f? (j). However, Eq. (5) is apt 
to "premature". So an enhanced Richardson-Lucy iterative 
algorithm [8] is applied in this paper. 

The nonlinear physical constraints, including upper limit 
and lower limit, can be expressed as [20] 


{ £0) > u() > FG) = u(): © 
FONG) « BG) > FO (9) = 8(). 


Obviously, each element of F is less than the upper limit. 
So the lower limit only needs to be considered. For y-ray 
spectra, background is regarded as the lower limit b(/). 

Equations (5) and (6) consist a complete representation of 
DDM. Referring to Ref. [8], exponential factor p was intro- 
duced to DDM. The principles of enhanced DDM can be ex- 
pressed as follows: 

Step 1: for n = 0, set the iterative initial value f(0) = 
[; Lee Ju 

Step 2: setting the number of iterations L and cycle number 
R; 

Step 3: initializing the cycle number r — 1; 

ips a according to Eq. (4) to calculate and seek the solu- 
tion f% 

228 5: introducing parameter p, and f?) (i) = [/ 02 ())]", 

€ (1,2),¢=0,1,---,N-1; 

E 6: introducing physical constraints; 

Step 7: if r = R, then stop calculating. Otherwise, r = 
r + 1, continue to do Step 4. 


fe» GH = FOU) ), 6 


050202-2 


HIGH-RESOLUTION BOOSTED... 


IIl. RESULTS AND EVALUATION 


The difficulty evaluation in separation of doublets of IEC 
standard [21] is given in Table 1. Referring to it, the Gaussian 
function was used to simulate nine full-energy peaks. Their 
positions, heights and areas are listed in Table 2. The syn- 
thetic spectrum is shown in Fig. 1(a). From Table 2, the 
biggest ratio in height is 10:1, and the minimum interval of 
peak positions is less than 1/3 FWHM. All cases except for 
area ratio 1:100 in Table 1 are considered. The spectral line 
was reconstructed by applying DDM with 1000 (L = 50, 
R = 20, p = 1.8), 5000 (L = 100, R = 50, p = 1.8), 
10000 (L = 200, R = 50, p = 1.8), and 50000 (L = 500, 
R = 100, p = 1.8) iterations, respectively. Then, peak areas 
were calculated by accumulating counts for “isolated” singlet 
of reconstructed spectrum. The result of 50000 iterations is 
shown in Fig. 1(b). 


TABLE 1. Difficulty evaluation in separation of overlapped peaks 
of IEC standard. 


S1/S2 >3/2 FWHM ~FWHM <1/3 FWHM 
1/1 Easy Medium Very difficult 
1/3 Medium Difficult Very difficult 
1/10 Medium Difficult Very difficult 
1/100 Very difficult Very difficult Very difficult 


TABLE 2. Difficulty evaluation in separation of overlapped peaks 
of IEC standard 


Peak No. Position (Ch.) Height Area 
1 30 500 6266 
2 50 5000 62667 
3 60 1700 21306 
4 80 1700 21306 
5 84 1700 21306 
6 110 400 5013 
7 114 4000 50133 
8 125 400 5013 
9 210 400 5013 


Also, width of the reconstructed peak decreased with in- 
creasing number of iterations. Theoretically, a complete spec- 
trum consisting of ó-function (single pulse) can be obtained 
after a large number of iterations. And corresponding counts 
of pulse are peak areas. But a bigger iteration means more 
time and memory consumed. Thus, from the perspective of 
efficiency, iteration is stopped when doublets are completely 
separated. The peak areas are not the corresponding counts 
of pulse, but the sum of corresponding counts of a few “iso- 
lated" points. The results show that DDM can effectively de- 
compose the peaks overlapped by just 1/3 FWHM. It should 
be emphasized that this is an extremely difficult task. 

Table 3 lists the peak positions and areas of the DDM- 
reconstructed spectrum. For area ratio of 1:1 and peak po- 
sition interval of about 1/3 FWHM, the error of peak position 
is 2 channels; while it is 1 channel for area ratio of 1:10 and 
peak position interval of about FWHM. The peak area error 
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Fig. 1. (Color online) synthetic spectrum (a) and reconstructed spec- 
trum using DDM algorithm after 50000 iterations (b). 


is less than 4% for all cases. This indicates that positions and 
areas can be precisely obtained with DDM. 


TABLE 3. Reconstructed peak positions, channel contents and their 
errors in the spectrum shown in Fig. 2 


Peak Position Position Channel Contents 
No. (Ch.) error (Ch.) contents error (%) 
1 30 0 6593 0.52 

2 50 0 64273 2.5 

3 60 0 21621 1.48 

4 80 0 20556 3.52 

5 82 2 20614 3.25 

6 110 0 5073 1.2 

7 113 1 49717 0.83 

8 125 0 5018 0.1 

9 210 0 5014 0.02 


IV. DISCUSSION 


Actually, a measured y-ray spectrum usually contains 
noise and background. Therefore, it is necessary to study the 
influence of noise and background on DDM. 
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Fig. 2. (Color online) A Gaussian peak and its DDM reconstructed spectrum. (a) and (b) peak with random noise and the DDM result; (c) 


and (d), peak with gauss noise and the DDM result. 


A. Noise 


A Gaussian peak is added with random (Fig. 2(a)) and 
gauss (Fig. 2(c)) noises, in amplitudes of 3% and 1% of the 
Gaussian peak height, respectively. The results are shown in 
Figs. 2(b) and 2(d), respectively. One sees that DDM has a 
strong inhibitory effect on noise. Because the least-squares, 
smoothing methods, and iterative can solve ill-conditioned 
problem, and the Richardson-Lucy iteration algorithm is 
adopted, noises can be intangibly inhibited in the process of 
iteration. If the noise is not serious, a spectrum can be recon- 
structed without de-noising. 


B. Background 


A spectrum is added with straight-line, oblique-line, and 
step backgrounds, respectively. The DDM-reconstructed re- 
sults are shown in Fig. 3. The spectra were not reconstructed 
correctly. This indicates that background affected the accu- 
racy of reconstructed spectra, but researches show that if the 
response matrix can be obtained under the same condition 
of background, or the background can be removed from the 
measured spectrum and response matrix by applying the same 
way, the lines can be reconstructed accurately. 


C. Establishment of response matrix with GEANT4 


Experiments have shown that the response function of a 
spectrometer system is strongly dependent on energy of the 
incident y-rays. Heretofore, there are no unanimous theo- 
ries and empirical expressions, but we could obtain the re- 
sponse function experimentally. Ideally, the response func- 
tion shall be calibrated by using standard radioisotope sources 
to have strong single peaks distributed uniformly through- 
out the whole energy range. However, it is impractical to 
obtain so many standard sources and the y-ray peaks can- 
not be in such a uniform distribution. Fortunately, Monte 
Carlo method can correctly simulate response function of y- 
ray spectra [22, 23]. Above all, distribution and energy of 
radioactive source can be artificially controlled. 

GEANT [24] can be used for neutron, photon, electron, or 
coupled neutron/photon/electron transport, proton, etc., and 
is capable of calculating eigenvalues for critical systems. For 
photons, the code accounts for incoherent and coherent scat- 
tering, possibility of fluorescent emission after photoelectric 
absorption, absorption in pair production with local emission 
of annihilation radiation, and bremsstrahlung. Above all, the 
photon energy regime is from 1 eV to 100 GeV for Rayleigh 
and Compton effects, down to the lowest binding energy for 
each element for photo-electric and ionization, and down to 
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TABLE 4. The scaling factor and net areas of the NO. 104 gamma-ray spectrometer 


Model Standard contents 


U (ug/g) Th (ug/g) 
U 191.38 5.95 
Th 8.31 371.28 
K 4.29 7.97 
Scaling factors Th (10^ $/cps) 13.316 
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Fig. 3. (Color online) Spectra with backgrounds in (a) straight line, 
(b) oblique-line and (c) steps, and the reconstructed spectra. 


Net peak areas (cps) 

K (%) U Th K 
0.25 25.73 0 3.34 
0.58 1.94 27.882 1.674 
4.48 0.505 0.746 11.325 

7.438 0.396 U (107 S/cps) K (%lcps) 


10eV for bremsstrahlung. GEANT4 was employed in the 
work. 


In order to coincide with experiment, properties of the de- 
tector must be revised before simulation. 7*!Am (0.06 MeV), 
?"Co (0.122 MeV), ?Na (0.51 MeV), !37Cs (0.661 MeV), and 
?4Na (1.38 MeV, 2.75 MeV) were selected. The crystal di- 
mensions were adjusted to keep spike and full-energy peak 
efficiency of the measured spectrum coincidence with the 
simulation results. Then, according to the data of energy cali- 
bration, corresponding energy per channel was chosen to sim- 
ulate response function. It was adjusted according to multi- 
channel analyzer. And the analog data of each single peak 
were normalized as response function of this energy. The 
response function needs to be converted to response matrix 
for its use in DDM: measurement spectrum is convolution of 
input data and response function. So referring to the convo- 
lution process, response function is reflexed, zero filled, and 
gradually shifted to get a series of two-dimensional vectors 
which constitute response matrix. 


V. CALIBRATION OF NAI GAMMA SPECTROMETER 
BASED ON DDM 


DDM was used to calibrate the field NaI gamma spectrom- 
eter. Generally, spectrum stripping was used. The princi- 
ples of calibration can be seen in Ref. [25]. IED-3000A Nal 
gamma spectrometer (104#) was used. The dimensions of 
Nal crystal were 75 mm x 75 mm, with an energy resolution 
of 10.3% (661 keV). The measuring time was 3600s. Sub- 
stance contents of quasi-saturation model were listed in Ta- 
ble 4, which were measured by Analysis and Test Center of 
China Geological Survey. 


The measured spectrum of hybrid model was shown in 
Fig. 4(a). SNIP [26] was applied to deduct background. 
DDM was used to reconstruct the spectrum. The result 
after 50000 iterations is shown in Fig. 4(b). Calibration 
coefficients and net peak areas of 1.46 MeV (K), 1.76 MeV 
(?4Bi of U series), and 2.62 MeV (PTI of Th series) are 
listed in Table 4. 


Table 5 shows the elemental analysis results of a y-ray 
analysis standard by applying DDM and spectrum stripping 
method (SSM), It can be seen that the DDM results are of 
smaller deviations from the given value of U, Th and K, being 
respectively —1.02%, —4.19% and —2.41%, while the SSM 
results deviated by —8.06%, 10.44% and 8.84%, respectively. 
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Fig. 4. (Color online) The measured spectrum of the U model (a) and the reconstructed spectrum (b). 


TABLE 5. Results of No. 104 gamma-ray spectra with DDM and 
spectrum stripping method (SSM) 


Elemental contents Deviations (96) 


Type of y-rays 


Standard DDM SSM DDM SSM 
U(.76MeV)(ug/g 61.60 61.06 56.72 —1.02 —8.06 
Th (2.62 MeV) (ug/g) 183.93 176.21 203.13 —4.19 10.44 
K (1.46 MeV) (%) 2.49 255 271 —2.41 8.84 

VI. CONCLUSION 


Due to response function and iterative algorithm, DDM can 
be rarely restricted by shape of y-ray spectra, radioactive in- 
tensity. Besides, DDM can significantly pseudo-improve en- 
ergy resolution, precisely search peaks, calculate areas, and 
improve the performance of gamma spectrometer without 
hardware cost. It is worth mentioning that response matrix 
is the key of DDM. However, in practical, it is tedious and 
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